// Copyright 2012 Dorival de Moraes Pedroso. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.

// gemlab wraps the C++ library
//  NOTE: __IMPORTANT__ this package works on 64bit machines only,
//                      since Go's int is converted to C's long
package gemlab

/*
#cgo CXXFLAGS: -O3 -I. -I/usr/include/eigen3 -I/usr/lib/openmpi/include -I/usr/lib/openmpi/include/openmpi -I/usr/include/metis -pthread
#cgo CFLAGS: -O3 -I.
#cgo LDFLAGS: -L. -lstdc++ -lm -L/usr/lib/openmpi/lib -lmpi_cxx -lmpi -ldl -lmetis -lparmetis
#include "connectgemlab.h"
#include <stdlib.h>
*/
import "C"

import (
    "errors"
    "fmt"
    "unsafe"
)

// Uregions contains all data for unstructured regions
type Uregions struct {
    Tags   []int       // region tags [nreg]
    Maxs   []float64   // max areas/volumes [nreg]
    SFtags []int       // segments or facets tags [nsegfac]
    Regs   [][]float64 // region indicators (points) [nreg][ndim]
    Holes  [][]float64 // holes [nhol][ndim]
    Points [][]float64 // points [npts][ndim]
    SegFac [][]int     // segments or facets [nsegfac][connectivity]
}

// Sregions contains all data for structured regions
type Sregions struct {
    Tags          []int       // tags of regions [nreg]
    Nxs, Nys, Nzs []int       // number of divisions [nreg]
    Axs, Ays, Azs []float64   // coefficients for divisions [nreg]
    Nlx, Nly, Nlz []int       // flags for nonlinear divisions [nreg]
    Points        [][]float64 // points [npts][ndim]
    Conn          [][]int     // region connectivities [nreg][connectivity]
    Btags         [][]int     // boundary tags [nreg][nsides]
}

// Lregions contains all data for lines (linear elements such as rods, beams, drains)
type Lregions struct {
    Rtags  []int       // tags of rods [nreg]
    Jtags  []int       // tags of joints [nreg]
    Quads  []int       // quadratic lines?
    Points [][]float64 // points [npts][ndim]
    Conn   [][]int     // region connectivities [nreg][connectivity]
    Btags  [][]int     // boundary tags [nreg][2]
}

// Vtags collects data to tag vertices using their x,y,z coordinates
type Vtags struct {
    Tags   []int
    Coords [][]float64
    Tols   []float64
}

// UmeshGen generates unstructured meshes
// NOTE: le is optional
func UmeshGen(fn string, quad bool, nparts int, pfull bool, globalmax float64, rg *Uregions, le *Lregions, vt *Vtags) (err error) {

    // regions
    rndim, nsef, nreg, nhol, tags, sftags, maxs, regs, hols, npts, pts, nconn, conn, rerr := getRegions(rg)
    if rerr != nil {
        err = rerr
        return
    }

    // lines
    lndim, nle, le_rtags, le_jtags, le_quads, le_npts, le_pts, le_nconn, le_conn, le_btags, lerr := getLines(le)
    if lerr != nil {
        err = lerr
        return
    }
    if nle > 0 {
        if lndim != rndim {
            return errors.New(fmt.Sprintf(_gemlab_err1, rndim, lndim))
        }
    }

    // filename
    cfn := C.CString(fn)
    defer C.free(unsafe.Pointer(cfn))

    // quadratic and full partition flags
    var cquad, cpfull int
    if quad  { cquad  = 1 }
    if pfull { cpfull = 1 }

    // vertex tags
    nvt, vt_tags, vt_coords, vt_tols, verr := serializeVtags("UmeshGen", rndim, vt)
    if verr != nil {
        err = verr
        return
    }

    // generate
    status := C.umeshgen(cfn, (C.long)(cquad), (C.long)(rndim), (C.long)(nparts), (C.long)(cpfull),
        (C.double)(globalmax),
        tags,  maxs, nsef, sftags,
        nreg,  regs,
        nhol,  hols,
        npts,  pts,
        nconn, conn,
        nle, le_rtags, le_jtags, le_quads, le_npts, le_pts, le_nconn, le_conn, le_btags,
        nvt, vt_tags, vt_coords, vt_tols,
    )
    if status != 0 {
        return errors.New(fmt.Sprintf(_gemlab_err2, status))
    }
    return
}

// SmeshGen generates structured meshes
// NOTE: le is optional
func SmeshGen(fn string, quad bool, nparts int, pfull bool, rg *Sregions, le *Lregions, vt *Vtags) (err error) {

    // blocks
    bndim, nreg, tags, nxs, nys, nzs, axs, ays, azs, nlx, nly, nlz, npts, pts, nconn, conn, btags, berr := getBlocks(rg)
    if berr != nil {
        err = berr
        return
    }

    // lines
    lndim, nle, le_rtags, le_jtags, le_quads, le_npts, le_pts, le_nconn, le_conn, le_btags, lerr := getLines(le)
    if lerr != nil {
        err = lerr
        return
    }
    if nle > 0 {
        if lndim != bndim {
            return errors.New(fmt.Sprintf(_gemlab_err3, bndim, lndim))
        }
    }

    // filename
    cfn := C.CString(fn)
    defer C.free(unsafe.Pointer(cfn))

    // quadratic and full partition flags
    var cquad, cpfull int
    if quad  { cquad  = 1 }
    if pfull { cpfull = 1 }

    // vertex tags
    nvt, vt_tags, vt_coords, vt_tols, verr := serializeVtags("SmeshGen", bndim, vt)
    if verr != nil {
        err = verr
        return
    }

    // generate
    status := C.smeshgen(cfn, (C.long)(cquad), (C.long)(bndim), (C.long)(nparts), (C.long)(cpfull),
        (C.long)(nreg), tags,
        nxs, nys, nzs,
        axs, ays, azs,
        nlx, nly, nlz,
        npts, pts, nconn, conn, btags,
        nle, le_rtags, le_jtags, le_quads, le_npts, le_pts, le_nconn, le_conn, le_btags,
        nvt, vt_tags, vt_coords, vt_tols,
    )
    if status != 0 {
        return errors.New(fmt.Sprintf(_gemlab_err4, status))
    }
    return
}

// auxiliary functions /////////////////////////////////////////////////////////////////////////////////////////////////

func serializeVtags(msg string, ndim int, vt *Vtags) (_nvt C.long, _tags *C.long, _coords, _tols *C.double, err error) {
    if vt == nil {
        return
    }
    nvt := len(vt.Tags)
    if len(vt.Coords) != nvt {
        err = errors.New(fmt.Sprintf(_gemlab_err21, msg, nvt, len(vt.Coords)))
        return
    }
    if len(vt.Tols) != nvt {
        err = errors.New(fmt.Sprintf(_gemlab_err22, msg, nvt, len(vt.Tols)))
        return
    }
    coords := make([]float64, nvt * ndim)
    for i := 0; i < nvt; i++ {
        if len(vt.Coords[i]) != ndim {
            err = errors.New(fmt.Sprintf(_gemlab_err23, msg, len(vt.Coords[i]), ndim))
            return
        }
        for j := 0; j < ndim; j++ {
            coords[i*ndim + j] = vt.Coords[i][j]
        }
    }
    _nvt = (C.long)(nvt)
    _tags, _coords, _tols = (*C.long)(unsafe.Pointer(&vt.Tags[0])), (*C.double)(unsafe.Pointer(&coords[0])), (*C.double)(unsafe.Pointer(&vt.Tols[0]))
    return
}

func serializePoints(msg string, npmin int, P [][]float64) (ndim int, _np C.long, _p *C.double, err error) {

    // check number of points
    npts := len(P)
    if npts < npmin {
        err = errors.New(fmt.Sprintf(_gemlab_err5, msg, npmin, npts))
        return
    }

    // check space dimension
    ndim = len(P[0])
    if ndim < 2 || ndim > 3 {
        err = errors.New(fmt.Sprintf(_gemlab_err6, msg, ndim))
        return
    }

    // points
    p := make([]float64, npts * ndim)
    for i := 0; i < npts; i++ {
        if len(P[i]) != ndim {
            err = errors.New(fmt.Sprintf(_gemlab_err7, msg, len(P[i]), ndim, i))
            return
        }
        for j := 0; j < ndim; j++ {
            p[i*ndim + j] = P[i][j]
        }
    }
    _np, _p = (C.long)(npts), (*C.double)(unsafe.Pointer(&p[0]))
    return
}

func serializeConn(msg string, npts int, nclim []int, C [][]int) (_nc, _c *C.long, err error) {

    // number of points connecting regions
    nreg  := len(C)
    nc    := make([]int, nreg)
    nctot := 0
    for i := 0; i < nreg; i++ {
        nc[i] = len(C[i])
        if nc[i] < nclim[0] || nc[i] > nclim[1] {
            err = errors.New(fmt.Sprintf(_gemlab_err8, msg, nclim[0], nclim[1], nc[i], i))
            return
        }
        nctot += nc[i]
    }

    // connectivities
    c := make([]int, nctot)
    for i := 0; i < nreg; i++ {
        for j := 0; j < nc[i]; j++ {
            if C[i][j] < 0 || C[i][j] > npts-1 {
                err = errors.New(fmt.Sprintf(_gemlab_err9, msg, 0, npts-1, C[i][j], i))
                return
            }
            c[i*nc[i] + j] = C[i][j]
        }
    }
    _nc, _c = (*C.long)(unsafe.Pointer(&nc[0])), (*C.long)(unsafe.Pointer(&c[0]))
    return
}

func serializeBtags(msg string, nreg, nbtag int, B [][]int) (_b *C.long, err error) {

    // boundary tags
    b := make([]int, nreg * nbtag)
    for i := 0; i < nreg; i++ {
        if len(B[i]) != nbtag {
            err = errors.New(fmt.Sprintf(_gemlab_err10, msg, nbtag, len(B[i]), i))
            return
        }
        for j := 0; j < nbtag; j++ {
            b[i*nbtag + j] = B[i][j]
        }
    }
    _b = (*C.long)(unsafe.Pointer(&b[0]))
    return
}

func getRegions(rg *Uregions) (ndim int, cnsef, cnreg, nhol C.long, tags, sftags *C.long, maxs, regs, hols *C.double, npts C.long, pts *C.double, nconn, conn *C.long, err error) {

    // number of regions
    if rg == nil {
        err = errors.New(_gemlab_err11)
        cnreg = 0
        return
    }
    nreg := len(rg.Tags)
    cnreg = (C.long)(nreg)
    if nreg < 1 {
        err = errors.New(_gemlab_err12)
        return
    }

    // check
    if len(rg.Maxs) != nreg || len(rg.Regs) != nreg {
        err = errors.New(fmt.Sprintf(_gemlab_err13, nreg))
        return
    }
    nsef := len(rg.SFtags)
    cnsef = (C.long)(nsef)
    if len(rg.SegFac) != nsef {
        err = errors.New(fmt.Sprintf(_gemlab_err14, nsef))
        return
    }

    // basic data
    tags, sftags, maxs = (*C.long)(unsafe.Pointer(&rg.Tags[0])), (*C.long)(unsafe.Pointer(&rg.SFtags[0])), (*C.double)(unsafe.Pointer(&rg.Maxs[0]))

    // serialize regions
    ndim, _, regs, err = serializePoints("getRegions (regions indicators)", 1, rg.Regs)
    if err != nil {
        return
    }

    // serialize holes
    var hndim int
    hndim, nhol, hols, err = serializePoints("getRegions (holes)", 1, rg.Holes)
    if err != nil {
        return
    }
    if hndim != ndim {
        err = errors.New(fmt.Sprintf(_gemlab_err15, hndim, ndim))
        return
    }

    // serialize points
    var pndim int
    pndim, npts, pts, err = serializePoints("getRegions (all points)", 3, rg.Points)
    if err != nil {
        return
    }
    if pndim != ndim {
        err = errors.New(fmt.Sprintf(_gemlab_err16, pndim, ndim))
        return
    }

    // serialize connectivities of segments/facets
    nclim := []int{2, 2}
    np    := len(rg.Points)
    if ndim > 2 {
        nclim = []int{2, np-1}
    }
    nconn, conn, err = serializeConn("getRegions (segments/facets)", np, nclim, rg.SegFac)
    if err != nil {
        return
    }
    return
}

func getBlocks(rg *Sregions) (ndim int, nblk C.long, tags *C.long, nxs, nys, nzs *C.long, axs, ays, azs *C.double, nlx, nly, nlz *C.long, npts C.long, pts *C.double, nconn, conn, btags *C.long, err error) {

    // number of blocks
    if rg == nil {
        err = errors.New(_gemlab_err17)
        nblk = 0
        return
    }
    nreg := len(rg.Tags)
    nblk  = (C.long)(nreg)
    if nreg < 1 {
        err = errors.New(_gemlab_err18)
        return
    }

    // check
    if len(rg.Nxs)  != nreg || len(rg.Nys)   != nreg || len(rg.Nzs) != nreg ||
       len(rg.Axs)  != nreg || len(rg.Ays)   != nreg || len(rg.Azs) != nreg ||
       len(rg.Nlx)  != nreg || len(rg.Nly)   != nreg || len(rg.Nlz) != nreg ||
       len(rg.Conn) != nreg || len(rg.Btags) != nreg {
           err = errors.New(fmt.Sprintf(_gemlab_err19))
           return
    }

    // basic data
    tags          = (*C.long)  (unsafe.Pointer(&rg.Tags[0]))
    nxs, nys, nzs = (*C.long)  (unsafe.Pointer(&rg.Nxs[0])), (*C.long)  (unsafe.Pointer(&rg.Nys[0])), (*C.long)  (unsafe.Pointer(&rg.Nzs[0]))
    axs, ays, azs = (*C.double)(unsafe.Pointer(&rg.Axs[0])), (*C.double)(unsafe.Pointer(&rg.Ays[0])), (*C.double)(unsafe.Pointer(&rg.Azs[0]))
    nlx, nly, nlz = (*C.long)  (unsafe.Pointer(&rg.Nlx[0])), (*C.long)  (unsafe.Pointer(&rg.Nly[0])), (*C.long)  (unsafe.Pointer(&rg.Nlz[0]))

    // serialize points
    ndim, npts, pts, err = serializePoints("getBlocks", 4, rg.Points)
    if err != nil {
        return
    }

    // serialize connectivities
    nclim := []int{4, 8}
    if ndim > 2 {
        nclim = []int{8, 20}
    }
    nconn, conn, err = serializeConn("getBlocks", len(rg.Points), nclim, rg.Conn)
    if err != nil {
        return
    }

    // serialize boundary tags
    btags, err = serializeBtags("getBlocks", nreg, 2*ndim, rg.Btags)
    return
}

func getLines(le *Lregions) (ndim int, nle C.long, rtags, jtags, quads *C.long, npts C.long, pts *C.double, nconn, conn, btags *C.long, err error) {

    // number of lines
    if le == nil {
        nle = 0
        return
    }
    nreg := len(le.Rtags)
    nle   = (C.long)(nreg)
    if nreg < 1 {
        return
    }

    // check
    if len(le.Jtags) != nreg || len(le.Quads) != nreg || len(le.Conn) != nreg || len(le.Btags) != nreg {
        err = errors.New(fmt.Sprintf(_gemlab_err20))
        return
    }

    // basic data
    rtags, jtags, quads = (*C.long)(unsafe.Pointer(&le.Rtags[0])), (*C.long)(unsafe.Pointer(&le.Jtags[0])), (*C.long)(unsafe.Pointer(&le.Quads[0]))

    // serialize points
    ndim, npts, pts, err = serializePoints("getLines", 2, le.Points)
    if err != nil {
        return
    }

    // serialize connectivities
    nconn, conn, err = serializeConn("getLines", len(le.Points), []int{2, 3}, le.Conn)
    if err != nil {
        return
    }

    // serialize boundary tags
    btags, err = serializeBtags("getLines", nreg, 2, le.Btags)
    return
}

// error messages
var (
    _gemlab_err1  = "gemlab.go: UmeshGen: regions and lines must have the same space dimension. rndim=%d != lndim=%d"
    _gemlab_err2  = "gemlab.go: UmeshGen: C++ code failed with status = %d"
    _gemlab_err3  = "gemlab.go: SmeshGen: blocks and lines must have the same space dimension. bndim=%d != lndim=%d"
    _gemlab_err4  = "gemlab.go: SmeshGen: C++ code failed with status = %d"
    _gemlab_err5  = "gemlab.go: %s: minimum number of points is %d. np=%d is incorrect"
    _gemlab_err6  = "gemlab.go: %s: space dimension must be either 2 or 3. ndim = %d is invalid"
    _gemlab_err7  = "gemlab.go: %s: all points must have the same space dimension. %d != %d for point %d"
    _gemlab_err8  = "gemlab.go: %s: number of points in list of connectivity must be in [%d, %d]. nc=%d for region %d is incorrect"
    _gemlab_err9  = "gemlab.go: %s: connectivity point must be in [%d, %d]. con=%d for region %d is incorrect"
    _gemlab_err10 = "gemlab.go: %s: number of boundary tags must be equal to %d. nb=%d for region %d is incorrect"
    _gemlab_err11 = "gemlab.go: getRegions: all regions with unstructured data must be provided"
    _gemlab_err12 = "gemlab.go: getRegions: all regions with unstructured data must be provided"
    _gemlab_err13 = "gemlab.go: getRegions: size of Maxs and Regs slices must be equal to each other and equal to nreg = %d"
    _gemlab_err14 = "gemlab.go: getRegions: size of SegFac must be equal to the size of SFtags and equal to nsef = %d"
    _gemlab_err15 = "gemlab.go: getRegions: space dimension of holes and regions must be equal to each other. ndim: %d != %d"
    _gemlab_err16 = "gemlab.go: getRegions: space dimension of points and regions must be equal to each other. ndim: %d != %d"
    _gemlab_err17 = "gemlab.go: getBlocks: all regions with blocks data must be provided"
    _gemlab_err18 = "gemlab.go: getBlocks: all regions with blocks data must be provided"
    _gemlab_err19 = "gemlab.go: getBlocks: size of N#s, A#s, Nl#, Conn, and Btags slices must be all equal to each other == nreg"
    _gemlab_err20 = "gemlab.go: getLines: number of Rtags, Jtags, Quads, Conn and Btags must be all equal to each other"
    _gemlab_err21 = "gemlab.go: %s: number of coordinates must be equal to the number of vertex tags. %d != %d"
    _gemlab_err22 = "gemlab.go: %s: number of tolerances must be equal to the number of vertex tags. %d != %d"
    _gemlab_err23 = "gemlab.go: %s: all vertex coordinates must have the same space dimension. %d != %d"
)
